home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / linpack / zgefa.f < prev    next >
Text File  |  1996-07-19  |  3KB  |  112 lines
  1.       subroutine zgefa(a,lda,n,ipvt,info)
  2.       integer lda,n,ipvt(1),info
  3.       complex*16 a(lda,1)
  4. c
  5. c     zgefa factors a complex*16 matrix by gaussian elimination.
  6. c
  7. c     zgefa is usually called by zgeco, but it can be called
  8. c     directly with a saving in time if  rcond  is not needed.
  9. c     (time for zgeco) = (1 + 9/n)*(time for zgefa) .
  10. c
  11. c     on entry
  12. c
  13. c        a       complex*16(lda, n)
  14. c                the matrix to be factored.
  15. c
  16. c        lda     integer
  17. c                the leading dimension of the array  a .
  18. c
  19. c        n       integer
  20. c                the order of the matrix  a .
  21. c
  22. c     on return
  23. c
  24. c        a       an upper triangular matrix and the multipliers
  25. c                which were used to obtain it.
  26. c                the factorization can be written  a = l*u  where
  27. c                l  is a product of permutation and unit lower
  28. c                triangular matrices and  u  is upper triangular.
  29. c
  30. c        ipvt    integer(n)
  31. c                an integer vector of pivot indices.
  32. c
  33. c        info    integer
  34. c                = 0  normal value.
  35. c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
  36. c                     condition for this subroutine, but it does
  37. c                     indicate that zgesl or zgedi will divide by zero
  38. c                     if called.  use  rcond  in zgeco for a reliable
  39. c                     indication of singularity.
  40. c
  41. c     linpack. this version dated 08/14/78 .
  42. c     cleve moler, university of new mexico, argonne national lab.
  43. c
  44. c     subroutines and functions
  45. c
  46. c     blas zaxpy,zscal,izamax
  47. c     fortran dabs
  48. c
  49. c     internal variables
  50. c
  51.       complex*16 t
  52.       integer izamax,j,k,kp1,l,nm1
  53. c
  54.       complex*16 zdum
  55.       double precision cabs1
  56.       double precision dreal,dimag
  57.       complex*16 zdumr,zdumi
  58.       dreal(zdumr) = zdumr
  59.       dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
  60.       cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
  61. c
  62. c     gaussian elimination with partial pivoting
  63. c
  64.       info = 0
  65.       nm1 = n - 1
  66.       if (nm1 .lt. 1) go to 70
  67.       do 60 k = 1, nm1
  68.          kp1 = k + 1
  69. c
  70. c        find l = pivot index
  71. c
  72.          l = izamax(n-k+1,a(k,k),1) + k - 1
  73.          ipvt(k) = l
  74. c
  75. c        zero pivot implies this column already triangularized
  76. c
  77.          if (cabs1(a(l,k)) .eq. 0.0d0) go to 40
  78. c
  79. c           interchange if necessary
  80. c
  81.             if (l .eq. k) go to 10
  82.                t = a(l,k)
  83.                a(l,k) = a(k,k)
  84.                a(k,k) = t
  85.    10       continue
  86. c
  87. c           compute multipliers
  88. c
  89.             t = -(1.0d0,0.0d0)/a(k,k)
  90.             call zscal(n-k,t,a(k+1,k),1)
  91. c
  92. c           row elimination with column indexing
  93. c
  94.             do 30 j = kp1, n
  95.                t = a(l,j)
  96.                if (l .eq. k) go to 20
  97.                   a(l,j) = a(k,j)
  98.                   a(k,j) = t
  99.    20          continue
  100.                call zaxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
  101.    30       continue
  102.          go to 50
  103.    40    continue
  104.             info = k
  105.    50    continue
  106.    60 continue
  107.    70 continue
  108.       ipvt(n) = n
  109.       if (cabs1(a(n,n)) .eq. 0.0d0) info = n
  110.       return
  111.       end
  112.